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SHAPE AWARE MATCHING OF IMPLICIT SURFACES BASED ON THIN SHELL 

ENERGIES 

JOSE A. IGLESIAS, MARTIN RUMPE, AND OTMAR SCHERZER 


Abstract. A shape sensitive, variational approach for the matching of surfaces considered as thin 
elastic shells is investigated. The elasticity functional to be minimized takes into account two different 
types of nonlinear energies: a membrane energy measuring the rate of tangential distortion when de¬ 
forming the reference shell into the template shell, and a bending energy measuring the bending under 
the deformation in terms of the change of the shape operators from the undeformed into the deformed 
configuration. The variational method applies to surfaces described as level sets. It is mathematically 
well-posed and an existence proof of an optimal matching deformation is given. The variational model is 
implemented using a finite element discretization combined with a narrow band approach on an efficient 
hierarchical grid structure. Eor the optimization a regularized nonlinear conjugate gradient scheme and 
a cascadic multilevel strategy are used. The features of the proposed approach are studied for synthetic 
test cases and a collection of geometry processing examples. 


1. Introduction 

We present a variational model for the matching of surfaces implicitly represented as level sets. 
The approach is inspired by the mathematical theory of nonlinear elasticity of thin shells. The model 
consists in an energy functional, which is to be minimized among deformations of a computational 
domain in which two given surfaces are embedded. A minimizer of this functional is a deformation 
that closely maps one (reference) surface onto the other (template) surface. As the underlying model 
we consider the reference surface as a thin elastic shell, i.e., a layer of an elastic material embedded in 
a volume of another several orders of magnitude softer isotropic elastic material. Subject to matching 
forces the volume is deformed in such a way that the thin shell is mapped onto the template surface. 
The functional reflects desired phenomena like resistance to compression and expansion of the surface, 
resistance to bending, and rotational invariance, while solely involving the deformation and the Jac¬ 
obian of the deformation. The model is formulated in terms of projected derivatives from the tangent 
space of the reference surface onto the expected tangent space of the template surface. Taking into 
account a suitable factorization of the natural pullback under a deformation of shape operators enables 
us to formulate a model with appropriate convexity properties. The actual surface matching constraint 
is handled through a penalty, allowing for efficient numerical computation. 

Through arguments of compensated compactness, we are able to show weak lower semicontinu¬ 
ity of the energy and consequently existence of minimizing deformations. We present a numerical 
approach based on a multilinear finite element ansatz for the deformation implemented on adaptive 
octree grids. The resulting discrete energy is minimized in a multiscale fashion applying a regularized 
gradient descent. 

In the conference article ll32l a preliminary version of this approach was presented. For the func¬ 
tional in that paper lower semicontinuity could not be ensured for either the membrane or bending 
energies. This lack of lower semicontinuity manifests itself in applications, where compression of the 
surface is expected, and leads to undesired oscillations in almost-minimizing deformations, which we 
explore in the present work through explicit examples and computations. Additionally, to increase the 
efficiency the computational meshes are in the present paper adapted to the surfaces. Consequently the 
number of degrees of freedom scales asymptotically almost like that of a surface problem. 
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The main pillar of our modelling is the use of polyconvex energy densities, first introduced in lO. 
Energies of this type allow for geometric consistency properties like rotation invariance and the ability 
to measure area and volume changes. The core insight of this theory is that integrands consisting of 
convex functions of subdeterminants of the Jacobian give rise to integral functionals that are weakly 
lower semicontinuous in suitable Sobolev spaces. Indeed, this can be seen as an instance of com¬ 
pensated compactness ll50l . A generic polyconvex isotropic energy density of the type used in this 
work is 

(1.1) + /3,||Cof Af + r(det A), 

for Cof A:= det A A“^ the cofactor matrix of A. Here, the coefficients and the function T are such 
that ( |1.1| ) attains a local minimum for A G SO(n), that is, for rigid motions. Such an example is 
provided below in p.8| ). Often in the modelling of nonlinear elasticity the condition 

(1.2) lim r(A) = +00 

det A^0+ 

is added, to reflect the non-interpenetration of matter ||4l. In our model we make use of densities both 
with and without this property. 

Related work. Linear elasticity has been extensively used in computer vision and in graphics. Prom¬ 
inent applications are image registration |[48l[38l[55l|33l[34l, optical flow extraction (SSl, and shape 
modeling ll2^ . Recently, theories of nonlinear elasticity have been applied in many computer vision 
and graphics problems such as mesh deformation ifT^ . shape averaging llSbl . registration of medical 
images na. The advantage of nonlinear models is that they allow for intuitive deformations when the 
displacements are large. 

In this paper, we present a model for nonlinear elastic matching of thin shells. A finite element 
method for the discretization of bending energies of biological membranes has been introduced in Q. 
Their approach uses quadratic isoparametric finite elements to approximate the interface on which the 
gradient flow of an elastic energy of Helfrich type is considered. The papers (91 [TOl discuss accurate 
convex relaxation of higher order variational problems on curves described as jump sets of functions of 
bounded variation. In particular, it enables the numerical treatment of elastic energies on such curves. 

One challenge in polyhedral surface processing is to provide consistent notions of curvatures and 
second fundamental forms, i.e., notions that converge (in an appropriate topology or in a measure the¬ 
oretic sense) to their smooth counterparts, given a smooth limit surface. One computationally popular 
model for discretizing the second fundamental form is Grinspun’s et al. discrete shells model (301 . An¬ 
other efficient, and robust method for nonlinear surface deformation and shape matching is PriMo (6l. 
This approach is based on replacing the triangles of a polyhedral surfaces by thin prisms. During a 
deformation, these prisms are required to stay rigid, while nonlinear elastic forces are acting between 
neighboring prisms to account for bending, twisting, and stretching of the surface. We refer to Botsch 
and Sorkine (71 for a discussion of pros and cons for various such methods. In comparison with 
methods based on polyhedral surfaces, level set approaches like ours are not dependent on specific 
triangulations of the shapes. 

The matching of surfaces with elastic energies has recently been studied in (bH. Their energy 
contains a membrane energy depending on the Cauchy-Green strain tensor and a bending-type energy 
comparing the mean curvatures on the surfaces. The matching problem is formulated in terms of a 
binary linear program in the product space of sets of surface patches. For computations, a relaxation 
approach is used. 

A different direction is the use of parametric approaches to reduce shape matching problems to 
the matching of functions on a fixed domain. For example, the methods presented in (621 and (60l 
are based on conformal maps from the unit disk. A more general variant using conformal maps on 
surfaces with arbitrary topology is presented in (42l. Within the family of parametric methods, a 
surface matching approach related to ours is presented in ll44l . where nonlinear elastic energies are 
used for matching parametrized surface patches. In comparison to all these methods, our level set 
approach is non-parametric and allows surfaces of any topology, which does not need to be fixed in 
advance. 
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In ll59l , face matching based on a matching of corresponding level set curves on the facial surfaces 
is investigated. To match pairs of curves an optimal deformation between them is computed using an 
elastic shape analysis of curves. Compared to our approach, this model does not take into account 
bending dissipation of the curves. 

A different direction in shape recognition and matching is exploiting the intrinsic geometry of the 
surfaces only, thereby producing isometry-invariant methods based on the first fundamental form, like 
those in [l2^[m . In comparison, bending is penalized in our model and we use all curvatures of the 
surfaces and their directions to be able to better match regions of edges and creases correctly. 

A method for matching and blending of curves represented by level sets has been presented in 
l(49l . Thereby, a level set evolution generates an interpolating family of curves, where the associated 
propagation speed of the level sets depends on differences of level set curvatures. In this class of ap¬ 
proaches, geometric evolution problems are formulated, whereas here we focus on variational models 
for matching deformations. Variational registration of implicit surfaces was also considered in BOl . 
but only through volume elasticity, in contrast to our shell terms. 

To summarize, the main novelty of our contribution is the combination of independence of mesh 
topologies arising from the use of level sets, penalization of tangential distortion in a rotationally- 
invariant framework, and awareness both of curvatures and curvature directions of the surfaces in the 
matching. We are not aware of any other methods possessing all of these features simultaneously. 

Our approach is inspired by the articles ll2Tl[^ in which surface PDE models are derived in terms 
of the signed distance function. Shape warping based on the framework of ll2ll has been discussed 
from a geometric perspective in l[T4ll . 

Outline. The paper is organized as follows. In Section 2, we review the required preliminaries about 
distance functions and formulate the geometric non-distortion and matching conditions that inspire 
our model. In Section 3, we present the different contributions to our energy. Section 4 is devoted 
to proving the existence of minimizing deformations under suitable Dirichlet and Neumann boundary 
conditions. Furthermore, the strong convergence of solutions for vanishing matching penalty parameter 
is discussed and counterexamples showing the lack of lower semicontinuity of related simpler models 
are given. In Section 5 a numerical strategy for minimizing the energy on adaptive octree grids is 
presented. Finally, Section 6 contains a range of numerical examples demonstrating the behaviour of 
solutions corresponding to our design criteria, and presents several potential applications. 

Some useful notation. For later usage and the purpose of reference let us collect some useful notation, 
mostly introduced in detail in later sections: 

• \B\ stands for the Febesgue measure of S C and dmmB = \x — y\ for its 

diameter. 

• Generic matrices are denoted by A, S, M, N. We use II for the identity matrix. The set of 
rotations is denoted by 0(n) and SO(n) is the set of orientation-preserving rotations. The set 
of all symmetric and positive definite matrices is SPD(n). 

• Components of vectors are denoted with subindices. For v G It’I denotes its Euclidean 

norm. The {n — 1)-dimensional sphere is For a matrix M, |M| is the Frobenius norm. 

• For two column vectors v^w ^ 0 u; is the tensor product of v and w, that is, the square 

matrix vw^. In particular, if |i(;| = 1 we have the identity {v (g) w)w = v. 

• P(e) = II — e (g) e is the projection onto vectors orthogonal to e G 

• Deformations on R'^ are denoted by 0, and deformations defined on a hypersurface M d R^ 
by if. The identity deformation is denoted by id. 

• C R'^ denotes the computational domain. Every relevant deformation 0 maps into 12 
has to contain all computationally relevant manifolds M. 12 has Fipschitz boundary, is open 
and bounded. 

• We use the notation di for partial derivatives, V for the gradient of a scalar function, V for the 
Jacobian matrix of a vector function and 22^ for the Hessian matrix of a scalar function. 

• A4i, Af 2 C 12 are (7^’^ compact hypersurfaces. The inside and outside components of 12 \ M-i 
are well defined by the Jordan-Brouwer separation theorem (EH, Chapter 2, Section 5). 
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The signed distance function to Af i, A ^2 is denoted by di, d 2 . The sign convention is that 
d^ is negative on the inside of A4i, so that di(x) = — dist(x, A4i) if x is in the inside compon¬ 
ent of A4i and di(x) = dist(x, Mi) otherwise, where the distance functions dist(-, A4i), 
i = 1,2 are the unique viscosity solutions of 1 — |V dist(-, A4i)\ = 0 and dist(-, Mi) = 0 on 
Mi. The normal fields to the offsets of Mi at a point x are denoted by n^(x) := Vdi(x). A 
superscript next to Mi (i = 1, 2), as in M^, denotes that we are talking about a level set of d^ 
with value different from zero, so that M^ := d^^(c). 

denotes the tangent space to at x. The outwards normal to is 

given by ni{x), and the set of points where d^ is not differentiable is denoted by sing d^. 

We use Si = V‘^di for the Hessian of d^, which coincides with an extended shape operator 
of Mi. 

• A, /i are the Lame coefficients of an isotropic material in linearized elasticity. 

• K^) is the space of continuous functions from the domain Q to the range [R^, the 
Holder spaces in which the k-th derivative is a-Holder continuous, including the Lipschitz case 
(T = 1. The range of the spaces is specified unless it is R. Sobolev spaces are denoted by 

and the closure of compactly supported smooth functions in them by Wq'^. 

• The letter C is reserved for a generic positive constant that may have different values in each 
appearance. Sequence indexing is usually denoted by a superscript k, and limits by an overline, 
e.g., (j)’^ -> 

2. Deformation and matching of level set hypersurfaces 

We are given two compact, connected embedded hypersurfaces Adi, Af 2 of class which are 
diffeomorphic to each other, and both of which are contained in a bounded Lipschitz domain C R'^. 
In this section we deal with the tangential distortion and the change of the shape operator under a 
deformation 0 : 12 ^ 

For any c G K, we denote the c-offsets to the hypersurface Mi by M^ := {x G 12 | di(x) = 
c}. Furthermore, we define the singularity set sing as the set of points where d^ is not twice 
differentiable. With the regularity of Mi that we have assumed, it is well known (e.g.. Theorem 1.1, 
Corollary 1.3 and Remark 1.4 of (431) that sing d^ has Lebesgue measure zero and dist(A4^, sing d,)> 
0. Furthermore, combining (211 Theorem 5.6] and (45l Proposition 4.6, 7.] we see that d^ G (7^(12 \ 
sing di). 

The gradient of the signed distance function Vdi(x) is the outward-pointing unit normal ni(x) to 

dlx) d.(x) d.(x) 

M^ ^ at a point x. The tangent space to M^ ^ at x, denoted by T^M^ , consists of all vectors 
orthogonal to ni(x). Then, the corresponding projection matrices onto the tangent spaces are defined 
by 

Pj(x) := P(ni(a;)) = 1 - ni{x) ® r\i{x). 

Note that Si{x) := X>^dj(x) = Dnj(a;)Pj(x) is the shape operator of the immersed hypersurface 
A/|A(x) ^ point X. In fact, from |ni(x)p = 1 we deduce by differentiation that {x)Si{x) = 0. 

This, together with the fact that rii (g) rii is the projection onto the normal of the hypersurface Si shows 
that 

Vi{x)Vni{x) = Vni{x). 

With our choice of signs for d^, the symmetric matrices Si are positive semidefinite for convex hy¬ 
persurfaces Mi. Further information on tangential calculus for level set functions may be found in 
Chapter 9 of (251. 

2.1. Tangential derivative and area and length distortion. First, let us assume that (j) exactly maps 
M\ onto Ml, for all c> 0. Then, = im Pi(®) and = 

imP 2 (^(a;)) and we define the tangential derivative induced by the deformation 0 as 

Vig4>{x) V 2 {(l){x))V(j){x)Vi{x ), 


( 2 . 1 ) 
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capturing the tangential variation of (p{x) on M .2 along tangential directions on Mi. In the variational 
model we consider below an energy term depending on Vig(j){x) will reflect the tangential distortion 
of the deformation in the context of a matching of the two hypersurfaces Mi and M 2 even though 
4>(Mi) does not necessarily equal M 2 - Indeed, in the case M 2 ^ (l>{Mi) the variation along a tangent 
direction on AIi is still projected via Vtg(j){x) onto the tangent space not onto the 

tangent space of the deformed hypersurface 0(AIi) (cf. Fig. [^. Therefore there may exist tangential 

directions v G TxMi^^^\ such that = 0 even though Vc/iv ^ 0. Thus Vtg(j){x) can only 

be considered a measure of tangential distortion if (/)(Mi) is sufficiently close to M2 in the sense of 
closeness of tangent bundles. 



Figure 1. A sketch of the tangential derivative in the non-exact matching case 

with (/)(Mi) ^ M 2 - 

For a general deformation ^ ^ R'^ the Cauchy-Green strain tensor describes (up to 

first order) the deformation in a frame invariant (with respect to rigid body motions) way. Since we 
are interested in the effect of such a deformation between two hypersurfaces, for a suitably extended 
tangential gradient Ag0 + n 2 o 0 (g) ni we define the extended tangential part of the Cauchy-Green 
strain tensor, measuring only tangential distortion: 

(2.2) (Ag0 + (n 2 O 0) (g) m)^ (Ag0 + (n 2 O 0) (g) m) = T>tg0^Ag0 + ni (g) ni . 

The term u. 2 {f{x)) ® ni(x) is used to complement directions that are removed by the projections in 
the definition of the tangential distortion Ag0 and can be seen to realize a nonlinear Kirchhoff-Love 
assumption ifTTl Page 336], which postulates that lines normal to the middle surface of a shell remain 
normal after the deformation, without stretching. 

Next, we investigate the area and length distortion due to the tangential derivative Ag0- For a given 
vector e G we denote by Q(e) any proper rotation such that Q{e)en = e, where denotes the n-th 
element of the canonical basis of R^. Note that this condition does not specify a unique Q(e). Then, 
for every B ^ satisfying w G ker B and imB C for some unit vectors v^w ^ we have 

(2.3) Q{v)^{B + v® w)Q{w) = Q{v)^BQ{w) + (g) = 

where B is the upper left (n — 1) x (n — 1) submatrix of BQ{w) . Obviously p.3| ) implies 
det{B-\-v®w) — det(S) , 

\B + V ® wf — tr((S -h T’ (g) wY' {B -h T’ (g) w)) = 1 + \B\^ . 

Hence, for f(Mi) = M 2 and v = n 2 {f{x)), w = ni{x) the area distortion under the hypersur¬ 
face matching deformation f at some position x is described by det(Ptg0(^) + (g) ni(x)), 

which equals the positive square root of the determinant of the above Cauchy-Green strain tensor 
Ptg0^Ag0 + ni (g) ni. The squared tangential length distortion (in the sense of summing all squared 
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distortions with respect to an orthogonal basis) is described by \Vig(/){x) + n 2 ( 0 (x)) 0 ni(x))|^ and 
equals the trace of the Cauchy-Green strain tensor. 


2.2. Bending and curvature mismatch. Now, we quantify the change of curvature directions and 
magnitudes under the deformation 0. Our approach is motivated by models describing bending of 
elastic shells, because in our application the hypersurfaces are considered as thin shells. 

In order to quantify the changes of curvature we first assume that = M 2 , and compute the 

difference of the pull back of the shape operator S2 on M2 onto Mi under the deformation (j) and the 
shape operator Si on Af 1 , which, for two arbitrary directions t’, G is given by 

S2{<^{x))V(^{x)v • V(j){x)w — Si{x)v • w — {V<^{x)'^S2{<^{x))V(^{x) — Si{x)) v • w . 

If w are tangent vectors in TxMi, this difference describes the relative shape operator. 

We define the extended relative shape operator 

(2.4) Sreiix) ■■= 'D(l){xf S 2 {(t){x))V(l){x) - <Si(a;). 

For n = 3 and when 0 is an isometric deformation between Mi and M2 (that is, V(j){x) is an 
orthogonal mapping on TxMi for all x G Mi), Sr el appears in physical models for thin elastic shells 
in the context of the F—limit of 3D hyperelasticity ll28l . Even though we do not necessarily expect 
our deformations to be tangentially isometric, we use this ansatz to compare curvatures of level sets in 
deformed and undeformed configuration, respectively. The following calculations shed some light on 
the properties of SreV- 

P^(d 2 o (/))(a;) = X>((X>(/>)'^(n 2 o ^))(a;) 

(2.5) ^ 

= V(t>{xf'S2{(t>{x))V(t>{x) + {n2{(l){x)))^v‘^(t>^{x). 

k=l 


The assumption that 1 ) = M2 can be rewritten as d 2 o (j){x) = 0 for x G Mi. Let us assume 
that in addition d 2 o 0 is a distance function (that is | V(d 2 o 0) | = 1 ), then d 2 o 0 is again a distance 
function, and since d 2 o 0 = 0 it follows that the left hand side of p.5| ) is the shape operator of the 
hypersurface Mi. The first term in the right hand side is the pullback of ^ 2 . 

Let us remark that the appearance of a second fundamental form is consistent with Koiter’s non¬ 
linear thin shell theory |[36l, (iTl Section 11.1]. Regardless of whether d 2 o 0 is a distance function or 
not, ( |2.5| ) implies that 


V(^{x)^S 2 {(i){x))V(p{x) - Si{x) 

( 2 . 6 ) 


y^(n 2 ((/>(a;)))fcP^(/>^(x) + I>^(d 2 o ^ - di)(a;) 

k=l 

n 

'^{n2{((>{x)))kV‘^(()'^{x), 
k = l 


in case d 2 o 0 = di. In the next section, we use the extended relative shape operator to derive a 
variational model for the mismatch of curvatures. 


3. Energy functional 

Given two hypersurfaces Mi and M 2 our ultimate goal is to describe best matching deformations 
(j), which map Mi onto AI 2 as the minimizer of a suitable energy. Thereby, different energy terms 
will reflect a set of matching conditions for a volumetric deformation 0 : D ^ and without a hard 
constraint 0 (A/fi) — M 2 '> 

• A membrane deformation energy E’mem penalizes the tangential distortion measured through 

'Digcj). 

• A bending energy E’bend penalizes bending as reflected by the relative shape operator. 

• A matching penalty E’match ensures a proper matching of the two hypersurfaces Mi onto M 2 
via a narrow band approach. 

• A volume energy E’voi enforces a regular deformation on the whole computational domain D. 
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Our approach is based on level sets. Hence, we replace the integration over a single hypersurface, 
i.e., A4i, for the first three energies by a weighted integration over a narrow band of width a with 
0 < (7 < dist(Afi, singdi). To this end we will make use of a cutoff function rjcr G C^(IR) with 
/[R = 1 and supp rja = [—cr, a]. Additionally, rfo- is assumed to be even and strictly decreasing 

in [0, +oo). 

In what follows we introduce the four energy contributions separately. 

3.1. Tangential distortion energy. Picking up the insight gained in Section [XT] we formulate the 
membrane energy in terms of the length and area change associated with the tangential distortion 

^tg0‘ 

(3.1) £’mem[</>] = 5 [ ria{di{x))W + n2(0(a;)) (g) ni(a;)) dx, 

Jn 

where VP is a nonnegative polyconvex energy density vanishing at SO(n). The weight 6 reflects the 
proper scaling of the tangential distortion energy in case of a thin shell model with shell thickness 6. 
The energy ( |3.1| ) vanishes only on deformations (j) whose Jacobian matrix V(j){x) maps 

isometrically onto for every point x G suppr/cr o di. In consequence, both tangential 

expansion and compression are penalized. 

Let US remark, that the extension Ag0 + ^2 (0(^) ) ® ni {x) of the tangential derivative V^gcj) defined 
in ( |2.1| ) with rank n — 1 can degenerate or be orientation-reversing depending on the local configuration 
of and AJ 2 at x (cf. Figure]^ for examples). 



Figure 2. Configurations in which for the (obviously isometric) identity we have 
det (Vtgl{x) + n 2 (llx) 0 ni(x)) =0 (left) and det (Vtgl{x) + n 2 (llx) (g) ni(x)) < 

0 (right) and thus the extended tangential derivative degenerates or reverses orienta¬ 
tion. 

Furthermore, the energy density W should not satisfy W(B) 00 for det B ^ 0. A straightfor¬ 

ward modification of the arguments of Ciarlet and Geymonat (HS), ifTbl Theorem 4.10-2) leads to a 
smooth integrand W which has isometries as local minimizers, with the correct invariance properties, 
and with a Hessian for S = II which matches the quadratic energy integrand of the Lame-Navier model 
of linearized elasticity. With given Lame coefficients A, /i > 0, we select the energy 

W{A) = 1\A\^ + ^^ (det . 

This density fits into the notation of (O), if we choose p = q = 2 and r{t) = ct‘^ -h de . 

3.2. Bending energy. Now, we discuss a variational formulation of the curvature matching condition 
V(j){xyS 2 {(l){x))V(j){x) — 5i(x), which is equivalent to a vanishing relative shape operator (cf. ( |2.4| )), 
where Si = = PiP^d^P^ for i = 1,2 are the shape operators on the hypersurfaces Mi and 

M 2 , respectively. At first sight, it appears natural to formulate a quadratic penalization and to define a 
bending energy 

Ebend[4’] = ^^ [ r]aidi{x))\V4>{x)'^S2{4>ix))V(l){x) - Si{x)\'^dx . 

Jn 
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The weight 6^ reflects the scaling of the bending energy for thin shells of thickness 6. However, 
this energy is in general not weakly lower semicontinuous. Indeed, consider a situation in which 
5i(x) = S 2 {(/){x)) = P(e^). A short computation shows that the corresponding density is not convex 
in its matrix variable along the rank-one segment joining 1 — |ei(8)ei and 1 — ^ei ® ei, which 
precludes lower semicontinuity (cf. also Example |4.8| below on the lack of rank-one convexity). In 
fact, this kind of density is closely related to the Saint Venant-Kirchhoff energy, whose quasiconvex 
envelope is computed in |[39]| . 

Thus, we are asking for an alternative lower semicontinuous energy functional which gives prefer¬ 
ence to deformations 0 for which o (l))V(j) is close to Si. We show that this can be achieved 

with the extended shape operators -h (g) for z = 1,2 and factorization. For 

proving this we make use of the following lemma. 

Lemma 3.1 (modified curvature matching condition). Assume that M, N are two symmetric, positive 
definite matrices satisfying 

M = PiMPi -h ni (g) ni and N = P 2 A^P 2 + n 2 (g) n 2 . 

Moreover, assume that A ^ satisfies 


^Pi=P2^, 

then the following statements are equivalent: 

(3.2) ^^P2iVP2^ = PiMPi 
and 

(3.3) A[M, N,A]:= P 2 iV^P 2 APiM“^Pi -h n 2 (g) ni G 0(n). 

Proof. By definition, the matrix A[M, W, A] is orthogonal if A[M, W, A]^A[M, A] = 1. Therefore, 
if p.3| ) holds, then 

II = (P 2 2P 2 APiAf 2P]^ -|- n2n^)^(P2 2 P 2 APiAf ^Pi + n2n^) 

= PiM“2P^^^P2W2P2P2W2P2APiM“ 2P^ + nin^n2nf 
= PiAT"^ A^P2(P2A^^P2)^P2AAr“^Pi + ninf . 

If we multiply this equation from left and right by PiAf sP^ and take into account that (P 2 A^ 2 P 2 )^ = 
P 2 A^P 2 and (PiAf 2 P ]^)2 = P^AfPi we see that this is equivalent to 

PiA^P2A^P2APi = PiATPi. 

Applying that APi = P 2 A we finally achieve at the equivalent condition 

A^P2A^P2A = PiATPi. 

The proof of the converse follows the same steps in opposite direction. □ 


If the assumptions of this lemma apply to Af = 5f^^(x), N = and A = V(j){x) with 

y = 0(x), then the curvature matching condition 

(3.4) (P</>(x))^P2(y)<Sr*(y)P2(y)P</>(x) = Pi(x)<Sr*(x)Pi(x) 


^^^{x)^S 2 ^^{(j){x))^V(j){x)) from 0(n 


and a lower semicontinuous energy func- 
would be a proper choice 


is equivalent to A{Sl^^{x),S 2 ^^{(j){x)),V(j){x)) G 0(n) 
tional penalizing deviations of A(5f 
for realizing curvature matching. Unfortunately, the positive definiteness assumption of Lemma ^A_ 
is not fulfilled if principal curvatures of Af 1 or M 2 are negative. Hence, we are replacing the exten¬ 
ded shape operator matrices by symmetric and positive definite curvature classification matrices 
Ci = C(S[^^), z = 1, 2, respectively. 

We have experimented with two different choices for C: 




SHAPE AWARE MATCHING OF IMPLICIT SURFACES BASED ON THIN SHELL ENERGIES 


9 


• A simple choice is + /il, where —fi is a strict lower bound of the principal 

curvatures. But in applications surfaces are frequently characterized by strong creases or rather 
sharp edges, leading to very large /i. As a consequence the relative difference of the eigenval¬ 
ues is significantly reduced when dealing with the resulting curvature classification matrices. 
Thus, the variational approach is less sensitive to different principal curvatures of the input 
hyper surfaces. 

• Another option is to use a truncation of the absolute value function for the eigenvalues of 
symmetric matrices. For a symmetric matrix B G with eigenvalues Ai,..., and a 
diagonalization B = Q^diag(Ai,... , \n)Q we use the classification operator 

C{B) = Q'^diagdAilr, • • •, \K\t)Q , 

where |A|r = max{| A|, r} for some r > 0. This approach properly represents the exact shape 
operator matching objective in case of principal curvatures of equal sign and absolute value 
larger than r. A disadvantage of this construction is that it is not able to force the deformation 
to correctly match curvature directions on the hypersurface with the same absolute value of the 
principal curvatures but with different signs. That is, locally a saddle point of the hypersurface 
may be mistaken for an elliptical point. However, this effect is usually compensated globally, 
and in applications the ansatz performs well, in particular in matching regions of edges and 
creases (see Sectionj^. 

Like for the membrane energy ( |3T| ), if V(j){x) is ensured to be orientation-preserving {deiVcj) > 0) 

and ni • (n 2 o 0) > 0 (cf. Figure |^, the curvature matching condition is equivalent to 

A{c{sr\x)),c{ST\m)),Bm) e son. 

Based on these considerations, a suitable choice for the bending energy is 

(3.5) £;bend[0] = / rjMi{x))W{K {C{Sl^\x)),C{St^\<i>{x))),V(l){x))) dx, 

Jn 

where W can be chosen as the same polyconvex density already used for F^mem- 


3.3. Mismatch penalty and volumetric regularization energies. So far, we have defined tangential 
membrane and bending energies which quantify the appropriateness of deformations 0 : fi ^ 
in a narrow band around the hypersurface A4i. In the derivation of these energies we assumed the 
constraint 0(A^i) = A42- However, such a constraint would be very hard to enforce numerically. 
Thus we use a weaker mismatch penalty instead: 


(3.6) 


-Ematch[0] = -[ iVcrO di)|d2 o 0 - di 

^ Ja 


2 


dx , 


where 1 /z/ is a penalization parameter. 

Moreover, we aim for a regular deformation on the whole computational domain Q which is 
globally injective. This, in particular, prevents from self-intersections of the deformed hypersurface 
i). To achieve this we introduce the following volume regularization term based on a polyconvex 
density W that enforces orientation preservation 


(3.7) 


Eyo\[^] = 


W (770, Cof 770, det Vcj)) dx 
-hoo 


if det 770(x) > 0 for a.e. x, 
otherwise. 


, where 


(3.8) IL(770,Cof770,det770) = 1770+ 00 Cof 7700 + 7 , (det 770)“", 

with p > n, q > n, s > {n — l)q/{q — n), and with > 0 ensuring that the density W attains 

a local minimum when 770^770 = 1. As mentioned in the introduction, such an energy is weakly 
lower semicontinuous in VF^’^(i7; R'^) when restricted to deformations whose Jacobian determinant is 
positive almost everywhere, and this condition is closed under weak convergence. 
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Choosing and using that because of its symmetries W can be expressed in 

terms of singular values ||20l Proposition 5.31], elementary but lengthy computations yield the station- 
arity condition at II 

(n + 1) {ap + {n - l)/3q) = ri^js, 

and that the corresponding Hessian is positive definite. For n = 3 an adequate example is then p = 
g = 4, s = 9, Op = 1, = 1, and 7 ^ = 4. For n — 2, one can use p — q — s — A, ap — Pq — 2 

and 75 = 3\/2. Notice that in this case, |P(/)| = |Cof P(/)|. 

3.4. Total energy. Sununing the above terms, our energy for shape-aware level set matching reads 
(3.9) Eu\(j)\ — -Ematchi^] ^mem[0] ^bend[0] 

where the different terms depend on the fixed input geometries Mi and M 2 through di and d 2 . 


4. Existence of optimal matching deformations 

First we prove the following weak continuity lemma, which is a generalization of the classical 
result given in |[50l Theorem 4.1]. Here the coefficients may depend on the deformed configuration. 

Lemma 4.1. Let ^ 0 G R'^) andp > n. Moreover, let Vi G x R^; i = 1, 2 

and we denote 

:= and Vi := </>(•)) , i = 1,2 . 

Then 

(4.1) det (P(vf )2?</)^P(v^) + vf ® v^) ^ det (P(v2)2?0P(vi) + V 2 ® vi) . 

_ 1 _ 

Moreover, for every symmetric positive definite Mi, i = 1, 2 with ^ G x R'^] and 

1 _ 

G (7^(12 X R'^] R'^^'^) and the corresponding compositions 

Mp-) and Mi := Mi{; ((>{■)) 


we have 


det A(Mf, M^,V<f)k) = det (P(v^)(M 2 ^) 5 P(v^)T>(/)^P(v^)(Mi^)- 5 P(v^) + ® v^) 

( 4 - 2 ) . _ ^ ^ _ 

-^det (P(v 2 )(M 2 ) 2 P(v 2 )X></)P(vi)(Mi) 2 P(vi) + V 2 ® vi) = det A(Mi, M 2 , . 

Proof. To prove ( |4.1| ) let Q € (fl). We show that 

:= [ Cdet(P(vf)P(/.'^P(v^)+v^( 8 )v^)dx ^ I := [ Cdet (P(v 2 )P(/)P(vi) + V 2 0 vi) dx. 

Jq Jq 

Moreover, we denote 


/ C det (P(v 2 )X>(/>^P(vi) + V 2 (8) vi) dx . 
Ja 

Using the inequality (cf. llTTl Theorem 4.7]) 

I det A-det 5| < C\A - B\max{\A\,\B\)^-^ 
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and Holder’s inequality it follows that 




< C [ \C\ P(v2)P0^P(Vi) - P(v2)P0^P(vi) + V 2 0 - V 2 0 Vi 

JQ 


max |P(vf + vf (g) I, \P{v 2 )V<p^P{vi) + V 2 ® viI 


n—1 


dx 


< c\ 


' Lp-n 


7-^ 


k\n-l 


+ 1 


P 

Ln-1 


> — V2 (8) Vi 


LP 


P(yi)V4>‘P{yf ) - P(v2)I70*P(vi) + vj 
< C||<||^^(||I7^.‘ing + l) 

■ jlP/llw (l|P(v5)lk»l|P(v?) -P{v,)|U» + l|P(vi)||i»||P(v‘) -P(V 2 )||l«) 

+ {||Vi - + ||V2 - V2 ||l‘«) 


Here, we have used that 


|P(v2)P0^P(vi) + V2 0 ViP <(\V(j)^\ + lY ^ <C(\V(j)^\^ + 1) • 

By the Rellich-Kondrakov embedding theorem (|[T1, Theorem 6.3 III) there exist subsequences of v^, 
i = 1, 2, which for simplicity of notation are again denoted by v^, i = 1,2, that converge uniformly to 
v^, i = 1,2, respectively. Taking into account the Lipschitz continuity estimate 

|P(e) - P(/)| = |(e - /) (8) e + / ® (e - /)| < 2^/n\e - f\ 

and that ^r^ —>■ Vj, i = 1,2 in L°° we obtain \I^ — I | —>■ 0 for /c —>■ oo. 

Next, we replace Vj, i = 1,2 in P by a piecewise constant approximation on a grid superimposed 
to the computational domain ft. Explicitly, we consider the finitely many non empty intersection 
cc| = S(z + [0,1]^) n of cubical cells with for 2 ; G and define 


fk 


Y] / Cdet (P(v2(2:^))P^^P(vi(2;^)) + V2{zs) 0 vi(2;^))da;, 

Jojf 


where zs is any point in n 6 c| if this set is nonempty. Using analogous estimates as above we obtain 


yk yk 

Is — I 


< 


ciia^j,(iii>/ii2g + i) 


(||P(vi,5)||loo||P(v2,5) -P(v2)||l°° + ||P(vi)||loo||P(vi^^) -P(vi)||j 


+ 


(||V2,, 


5 — V 2 IIL 00 + ||vi^5 — Vlllic 


0 . ’ 


Sluji = 


where vi <5 and \- 2 ,s are piecewise constant functions in with = vi( 2 ;( 5 ) and V 2 , 

\' 2 {zs), respectively. 

Using the uniform continuity of V 2 and vi on ft we obtain that 

ically increasing continuous function /3 : ^ K with /3(0) = 0. In particular the convergence is 

uniform with respect to k. The same argument applies for the difference of / and 


yk yk 


< f3{6) for a monoton- 


I 5 := / Cdet (P(v 2 ( 2 : 5 ))X><?^P(vi( 2 : 5 )) + V 2 ( 2 : 5 ) (g) Vi( 2 : 5 ))da: 

and we get |/^ — l| < C/3{S). Using ( | 2 . 3 | ) it follows that 

Q(v 2 ( 2 ;^))^(p(v 2 ( 2 ;a))AP(vi(z 5 )) + V 2 ( 2 ; 5 ) (g) Vl( 2 : 5 )j( 3 (vi( 2 ; 5 )) = 


' A 

0 ' 

0 

1 


Thus det(P(v 2 ( 2 ;( 5 ))AP(vi( 2 ;( 5 )) + V 2 ( 2 ;( 5 ) ( 8 )vi( 2 ;( 5 )) = det( 74 ) represents an (n — 1) x (n — 1) minor of 
the linear mapping corresponding to the matrix A with respect to different orthogonal basis in preimage 
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space (associated with P(vi( 2 :( 5 )) and vi( 2 ;( 5 )) and the image space (associated with P(v 2 ( 2 :( 5 )) and 
\' 2 {zs)). Indeed, denoting Qi := Q{'Vi{zs)) we have 

/ C(a:) det {P{v2{zs))'D(f)’^{x)P{vi{zs)) + V2{zs) ® vi(2;^))da: 

Ju^i 

= [ C{x)det{Q2{Pi-V2{zs))'D(f)’^{x)P{vi{zs))+ v2{zs)0vi{zs))Qi)dx 

= [ C{x) det {Q2Pi'V2{zs))Q2Q2'^(t>’"{x)QiQjP{vi{zs))Qi + e„ ® e„)da; 

Ju^l 

= I C(a^) det (p(en )<32 P 0 ^(a;)QiP(en) + Cn ® en)dx 

JujI 

= [ C(<3iy)det(P(e„)X)((3^o</)^oQi)(y)P(e„) + e„®e„)dy 

= [ C(<3l2/)Cofnn(^^(<32 0(/>^oQi)(y))dy, 

J OTcuf 


where we have used the orthogonal change of variables y = Qjx and Cofdenotes the minor ob¬ 
tained by erasing the last column and the last row. This change of orthogonal coordinates is fixed on 
each cell Since for each S the domain fl is covered by finitely many cells cc|, using the above com¬ 
putation and standard weak continuity results [|20l Theorem 8.20] for determinants of minors of the Jac¬ 
obian we obtain that ^ Is for k ^ oo. Finally, for given e we first choose S small enough to ensure 
that \ls — /| + |/| — I^\ < |. Then we choose A; large enough to ensure that |/^ — — Is\ < §• 

This proves that a subsequence of converges to / for A: ^ oo. Since the limit does not depend on 
the subsequence, we finally obtain weak convergence for the whole sequence. 

To prove ( |4.2| ), consider the three sequences of matrix functions 

(4.3) P(v^)(M|)ip(v^)+v^®v^, P(v^)X)</)^P(vJ^)+v^®vJ^ and P(vJ^)(Mf)"^P(vJ^)+vJ^®v^ 


The determinant of the second expression above converges weakly as A; ^ oo by the first part of the 
lemma, while the determinants of the first and third can be assumed to converge uniformly. Moreover, 
the matrices in ( |4.3| ) have the block structure shown in ( |2.3| ), so multiplying the three together and 
taking into account that P is a projection (depending on the argument) recovers the matrix 


P(v2")(M|)2P(v^)D0''P(vf)(M 


'=^-5p(v' 


+ ^2 


'^1 


appearing in the statement. Multiplicativity of the determinant and the fact that a product of strongly 
converging and one weakly converging sequence converges weakly then finishes the proof. □ 


We are now in a position to prove existence of a minimizing deformation for the hypersurface 
matching energy E’ in a suitable set of admissible deformations. Of particular difficulty is that deriv¬ 
atives of d 2 are not defined in the whole of Q and that in the functional these derivatives are evaluated 
at deformed positions. We handle this by ensuring that the involved deformations are such that terms 
involving these derivatives are not evaluated near the singularities. We obtain the following theorem: 

Theorem 4.2 (Existence of minimizing deformations). Let A4i, A42 be (7^’^ compact embedded hy¬ 
persurfaces in such that a diffeomorphism (^ : A4i ^ A42 exists between them. 

Assume further that 

(4.4) 0 < (7 < min(dist(A4i, singdi), dist(A42, singd 2 )), 

where singd^ is the set of points where d^ is not differentiable, and that C : ^ SPD(n) is 

continuous. Then there exists a constant 0 < uq := uo{Q^Mi^M 2 ^o-^p^ap) such that for 0 < ly < 
the functional has at least one minimizer 0 among deformations in the space VFq’^(12 ; [R’^) -h id. 
Moreover, 0 is a homeomorphism of 12 into 12, and (/)~^ G VF^’^(12; K^), where 9 is given hy 9 = 
q{l + s)/{q + s). 
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Proof. We proceed in several steps. 

Step 1: Coercivity. First we point out the coercivity enjoyed by our functional. Using the Poincare and 
Morrey inequalities ( BTIl . Theorem 12.30 and 11.34), and the Dirichlet boundary conditions we have 

(4.5) 11011(^0,a < C(1 + ||2^0 ||lp(O)) < C{1 + £’j,[(/)]p), 

for any f G + id and o = 1 — n/p. 

Step 2: Lower semicontinuity along sequences of constrained deformations. For the remainder of the 
proof, a deformation f G K^) + id, p > n is termed p-admissible for p > 0, if 

• Eyol[(p] < +00, 

• det V(j){x) > 0 for a.e. x G U, and 

• for all X G supp {rjcr o di) and every y G sing(d 2 ), we have |0(x) — y\'> p. 

Notice that since p > n, f has a unique continuous representative, so the third property is well defined. 
First, notice that with the assumption ( |4.4| ) we have 

(4.6) snpp{r]c 7 o d^) = {|di| < a} C U \ sing(di) , z = 1, 2. 

Let be a sequence of p-admissible deformations with Eyo\[f^] < C. By ( |4.5| ) and using the Banach- 
Alaoglu and Rellich-Kondrakov theorems, a subsequence (again denoted by (0^)) converges to a de¬ 
formation 0, both in the VF^’^-weak and uniform topologies. 

Now, we have (|[20l Theorem 8.20]) 

(4.7) (detr>/,Cof2?</)^) ^ {detV(f>,CofV(f>) inLu(fi) x . 

Additionally, since ( |4.7| ) holds and because Ejy[fk] is bounded, J^{detVfk)~^dx is bounded by 
the definition of W and det ^ 0 a.e. Together with we have 

(4.8) detV(l){x) > 0 a.e., 

SO that f is again p-admissible. 

Notice also that by a.e. positivity of the determinants, ( |4.7| ) and a standard lower semicontinuity 
result for convex integrands (see e.g. |[20l Theorem 3.23]) implies 

Eyoi[4>] < liminf E'voil^^], 

k^oo 

and uniform convergence of c/)^ immediately leads to 

^match[0] — lilD. F/matchi^ ]• 
k^oo 

We claim that under the assumptions of this theorem, we also have that 


(4.9) 

-E'mem [ 0 ] ^ HlH inf -£/mem 


k^oo 

and 


(4.10) 

Ehmd [f] < lim inf E’bendl^^] 
k^oo 


To see this, notice that 0^, f being p-admissible ensures that the normal vectors satisfy 

ni, n 2 o 0^, n 2 o 0 G 67^({|di| < a}; R'^). 

Consequently, the first part of Lemma |4 .1 1 (with Vi = n^) implies 

^{|di|<(T} det (Ptg0^+(n2 o 0 ^) (g) ni)) 

(4.11) ) \ ^ 2 V 

X{|di|<<T}(Ag0,det (r>tg</>+(n 2 o </)) gni)j in (L^’(f2))” x Ln{VL), 

with X{|dp<cr} denoting the indicator function. Combining ( |4.11| ) with the polyconvexity of VF, defin¬ 
ing the Action F^mem, both introduced in ( |3.1| ) we find the assertion ( |4.9| ). 

Furthermore, by our assumptions on Mi (see section]^, we have that 

e C0(n; R-X-). 
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Since C produces uniformly positive matrices, we have G We can 

then use a continuity result for square roots of nonnegative definite matrix-valued functions defined on 
Q |[T5l Theorem 1.1] to see that 


The second part of Lemma [4T] implies the weak convergence 

^{|di|<a}(A(C(5f*(<^'=)),C(5r‘(<^'=)),P<^'=),detA(C(5r*),C(5r*^ 


from which ( |4.10| ) follows by using the polyconvexity of W. 

Step 3: Existence ofminimizers restricted to admissible deformations. Since we have already seen that 
the set of p-admissible deformations is weakly closed and weakly compact, and that every term of E 
is weakly lower semicontinuous on this set, we just need to check that for all fixed u > 0, the set of 
p-admissible deformations, with adequate p, is not empty. 

For some given a satisfying dist(A42, sing d 2 ) — a > 0 let p satisfy 


(4.12) 


0 < p < dist(A42, singd 2 ) — a . 


We construct a deformation p, which is p-admissible and satisfies Ejy[(p] < oo. By assumption, there 
exists a diffeomorphism p : Mi ^ M 2 - Thus, we construct an extension of this diffeomorphism to 
{|di| < a} along the normal directions using 


(4.13) (p{x + sni{x)) (p{x) + 5n2(v:^(x)), for x G Afi, —a < s < a. 


We can then extend (f to the inside and outside components 12^, fio of 12 \ {|di| < cr} by solving 
the minimization problems for Eyo\ with Dirichlet boundary conditions given by ( |4.13| ) on dl2i and 
dQo \ dft, and by (p{x) = x on dl2. For the resulting p we have 


^match[^] — 0, £’vol[^] DC, F/memi^] DC, F/bend[^] ^ DC, 


where the first two statements follow by construction, and the last two by virtue of (p being a dif¬ 
feomorphism and the choice of a. Moreover, we note that since p has finite energy and the growth 
conditions assumed for W (see ( |3.7| )), the condition det V(f{x) > 0 for a.e. x is also satisfied JH. 

Step 4: A priori estimate to remove the constraint. Next, we show that for any p satisfying ( |4.12| ) there 
exists a parameter z/q > 0 such that for all 0 < z/ < z/q the constrained minimizers of Ei, subject to 
solves the unconstrained optimization problem, consisting in minimizing Ey on VFq- h id. 

To this end, we verify that every f that satisfies 


(4.14) 


Eiy[f] < E^[ 


is p-admissible. It is immediate from ( |4.14| ) that E’voK^) < +dc, and from the definition of W in p.7| ) 
it follows with the same arguments as in ( |4.8| ) that det 0 > 0 a.e. 

We prove now that for all deformations f satisfying ( |4.14| ) also satisfy 


(4.15) ||d2 0(/)||^oo({|di|<a}) < dist(Af2,singd2) - p . 

This is sufficient because from ( |4.15| ) it follows for all x satisfying |di (x) | < a by the triangle inequal¬ 
ity that 

p < dist(Af 2 , singd 2 ) - ||d 2 o (/)||Loo({|di|<a}) 

= dist(A42, singd 2 ) — dist((/)(x), A42) 

< dist((/)(x),singd2). 


which is the third property of a p-admissible deformation f. 

To prove (|4.15|) we use the triangle inequality and estimate 


||d2 o 0||Loo({|di|<a}) < cr -h ||d2 of- di||^oo({|di|<(j}). 


( 4 . 16 ) 
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By the monotonicity of rja- and the fact that the signed distance functions are Lipschitz continuous 
with constant 1 we have, for each a G (0, cr) that 

||d 2 o(j)- di||j:,oo({|di|<(T}) 

(4.17) < (l + ||0||co.“({cr-a-<|di|<(T}))+ 

Estimates ( |4.5| ) and ( |4.14| ) imply in turn 

ll0llcO-“({cr-5-<|dil<(T}) < < C{1 + Ei,[(p]p). 

Finally, combining ( |4.16| ), ( |4.17| ), and ( |4.18| ) we obtain 

||d 2 o 0||Loo({|di|<a}) 


11??^ o di(d 2 O 0 - di) ||L^({|di|<a-,7}) 
ricr{a - a) 


(4.19) <a + 


^1 + (7(1 + Ei,\(p\p^ 


0 -“ + 


-\\r]a o di(d2 o(j)- di)||ioo({|di|<^-a})- 


Voio- - d) 

Now we can apply Ehriing’s lemma l|54l Theorem 7.30] for the embeddings CC L°°{Q) C 

I/^(n) to control the last term in ( |4.19| ). Taking into account the Poincare inequality and Dirichlet 
boundary conditions, we obtain for any e > 0 a constant (7(e) > 0 such that 

W'qa O di(d2 0(j)- di)||ioo({|di|<a-a}) < her ° di(d2 O 0 - di)||j;^c«(Q) 

(4.20) < (7(e)||r?o- odi(d 2 o0 - di)||j;^ 2 (f 2 )+e C(||V(r?o- o di(d 2 o 0 - di))||ip(Q) + l) . 

Now, for the first term in the right hand side of ( |4.20| ) we can estimate 

(4.21) ||??<T0di(d2 0</)-di)||i2(j2) =I^^£^match[0]^ < E^[Lp\^ . 

For the second term, denoting diamfl = sup^ \x — y\, 

||V(77^ O dl(d2 0(j)- dl))||j;^p(Q) 

< \\V{r]a o di)(d 2 o(P- di)||i:,p(Q) + \\{ria o di)V(d 2 o ^ - di)||ip(Q) + 1 


< Cvp I ||d 2 o 


< 


dl||^L(j2)-®match[^]P j + C'(||27 (/>||/,p(q) + ||Vdi||iP(a) + l) 


(7i/p (^(||^||co,a(f 2 ) + 2diamfl) p E^[ip]p^ +c(^E^[ip]p + 1 j 


(4.22) 


< Cup {{1 + E, 


1 p-2 

Ip) p E^ 


Y^)+c(e,W + i), 


where we have applied the product rule, the definition of E’match, Ver ^ rja- < C, that |Vd^| = 1 
a.e., i = 1,2, the chain rule, and ( |4.14| ). The use of the chain rule is justified by ll46l Theorem 2.2], 
since d 2 has Lipschitz constant 1. 

Together, ( |4.20| ), ( |4.21D , and ( |4.22| ) imply 

\\r]a O di(d 2 O ^ - di)||Lc^({|di|«T-<T}) 

(4.23) ^ UP ^(7(6)12^ p E,^ \0\ 2 + e (7(1 + E^, [<^] p) p E^, [<^] p ^ + e (7 [<^] p + 1^ . 

In light of ( |4.19| ) and ( |4.23| ), and since Ej^[lp\ is independent of z/, we can now choose first a, then e 
and finally u small enough to obtain 

||d 2 o (/)||Loo({|di|<( 7 }) < cr + (dist(A42,singd2) - a - p) < dist(A42, singd 2 ) - p . 

Step 5: Injectivity. The injectivity and regularity of the inverse follow by the growth conditions satisfied 
by £Voi and classical results of Ball ||4l Theorems 2 and 3]. Note that Theorem 3 in [4| is stated in the 
mechanical application context in dimension n = 3, but it holds also in following the same proof 
and using the condition s > {n — l)q/{q — n). □ 


We have particularized the statement of Theorem |4.2| to the case of Dirichlet boundary conditions 
to ensure global invertibility. In fact, we also have existence of minimizing deformations for the case 
of Neumann boundary conditions. 
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Corollary 4.3 (Natural boundary conditions). Under the assumptions of Theorem |4.2| above, there 
exists a constant 

0 < I/JV = Ml, M 2 , CT,P, ap) 

such that for 0 < 1 / < i/jv, the functional E^, possesses at least one minimizer among deformations in 
the space R”). 


Proof. The proof follows the same arguments used for Theorem|4.2 


so we only point out the necessary 


modifications. We need a replacement for the coercivity estimate ( |4.5| ) and claim 

(4.24) ||0||vi/Tp(ri) ^ C{1 + + ||^0||LP(r2)) ^ C{1 + z/2£’^[0]2 + Ey[(j)\p). 

To verify this let us consider uo := {|di| < 6r/2}. An adequate Poincare inequality (see e.g. flTl 
Theorem 12.23]) implies that 

and we estimate the second term in the right hand side by 


/ 

J UJ 


(j)dx 


(/)dx< \(/)\dx< \d 2 0 (/)\dx 

J uo J uo Juo 

- di| dx + 

LU I 2 (z/£/pQatch [f] ) ^ “1“ 


+ \ijj\ sup \x\ 
xeM2 


< / 

d2 0 (j) 

J uo 



/(T\ -1 

b 

VI 

( 2 ) 


/ |di| dx + \(oo\ sup \x\ 
Juo xeA42 


f uo 

1 . _ . 1 


|di| dx + \(jj\ sup |x|, 
xeM2 


where Holder’s inequality has been used to compare and Lf norms. Therefore, ( |4.24| ) follows. 

The proof of the estimate for ||d 2 o (/)\\L°^(\\di\<a}) (lo ensure that deformations stay away from the 
singularities of d 2 ) is still valid with minor modifications, since u appears in ( |4.24| ) multiplicatively. 

□ 

We conclude this section with the following proposition, which explores the penalization limit in 
which the parameter o tends to zero. 

Proposition 4.4. Let {Tyk}ke\i^^ be a sequence of penalty matching parameters such that Ok ^ 0 SiS 
k ^ 00 , and be solutions of the Dirichlet minimization problem for Ei^j^. Then, up to a choice of 
subsequence, the converge strongly in to a minimizer of 

-E^mem £^bend ^'vol 

in K^) + II under the constraint 0(A4f) = for all c G (— cr, a). 

Proof First, notice that the energy E may be written as 

(4.25) EM= - f Vaodi\d2 0(t>-dif + ap\Vct>\P 
^ Ja 

+H ^ det Vif), CofX></), Vigcf), det(I>tg</> + n 2 o 0 (g) m), 

A(C(<Sf 0 o </>)), det (A(C(<SrO, C( 5 r* o </>)), P<^)) dx, 

where : [R+ x x x K x [R^x'^ x K ^ [R+ is smooth and convex. 

Denote by (p the extension of a diffeomorphism between Mi and M 2 used in the proof of Theorem 
Since S’matchl^] = 0, we have that Ej^^ [f^] < Ei[0\. By the coercivity estimate ( |4.5| ) the are 


4.2 


then bounded in and we may extract a (not relabelled) subsequence converging uniformly and 
weakly in to some limit f. Since {E^^^ [f^]} is bounded and Ok 0, the uniform convergence of 
implies that 


(4.26) [ 77^(di)|d2 

Jq 




di P dx 


k^OG 


[ Va{di)\d2 

Jq 


o 0 — dip dx = 0. 
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In consequence, C Since is the uniform limit of the maps which are sur¬ 
jective onto A ^2 *^1 is compact, we conclude that (^(A4f) = for all c G (—a, a). Therefore, 

(j) is admissible for all and [0^] < Ei[(j)]. Combined with lower semicontinuity and ( |4.25| ), the 
above implies 

(4.27) [ ap\V^^\P + H{det{V^^),...)dx -^ [ ap\V^\P + H{det{V^),...) dx. 

Jn Jn 

From this identity, the fact that H is convex and differentiable, and Vcj)^ Vcj) in it follows that 


0 = limsup ( / Op (^\V(j)^\^ — \V(j)\^\ + H{dei{V(j)^)^ • • •) ~ H{dei{V(j))^...) dx 
k^oo \Jq ^ ' 

> limsup ( f ap (\V(f>'^\P - \V(f>\p) + VH{det{V(f>),...) ■ (det(2?</)'=) - det{V(f>), 

k^oo \Jq ^ ^ 

= limsup f ap\'D(^^\'^dx — f ap\V(^\'^ dx . 
k^oo Jq Jq 

Together with the weak lower semicontinuity of the L^-norm, the above shows that 



f ap\V(p\'^ dx = lim f ap\V(^^\'^dx . 

Jn Jn 

Because LP(fl) has the Radon-Riesz property (||47l 2.5.26]), weak convergence and convergence of 
the norm guarantee strong convergence. Since was assumed to converge uniformly, we have also 
(/)k ^ (pin LP, and this shows that ^ 0 in R'^). 

That 0 is a minimizer of the constrained problem follows directly ([[8l, Theorem 1.21) from the 
fact that the Ejyj^ are an equicoercive family of functionals, T-converging in the weak topology of 
Indeed, equicoercivity follows easily from the above, while T-convergence is implied by the 
fact that Ejyj^ is an increasing sequence (|[8]|, Remark 1.40), because z/p ^ 0 appears as a denominator 
in F/match • n 


Remark 4.5. By the coercivity estimate ( |4.24| ) of Corollary |4. 3 [ an entirely analogous result holds for 
minimizers with Neumann boundary conditions. 


Remark 4.6. Contrary to what might be expected, the limit problem we have obtained is not a surface 
problem, since all the level sets are still coupled through the volume energy E’voi- The line of reasoning 
above depends heavily on the fact that the coefficients of the volume term are held fixed, since the 
equicoercivity and uniform strict quasiconvexity (in the language of f25i ) both require the presence of 
\\'^^\\^LP{n) functional. 

4.1. Oscillations and lack of rank-one convexity for the naive approach. To model the tangential 
distortion energy we have considered a frame indifferent energy density with the argument + 
(n 2 o 0) (g) ni. Let us now consider the case n — 2 and a simpler version of the membrane energy 
( [T6l ), where we use as an argument of the energy density directly the tangential Cauchy-Green strain 
tensor (cf ( |2.2| )) {T>tg(p{x))'^{'Dtg<p{x)) + ni(x) (g) ni(x), and define the membrane energy 

(4.28) ^mem[</>] := J^ricr{di{x))w(^{'Dtg(l){x))^Vtg(t){x) + ni(a;) ® ni(a:)j dx, 

with Vtgfp := 77(/)Pi defined as the tangential part of the derivative along and W : ^ 

R a frame indifferent energy density that has a strict minimum at SO(2). In fact, this energy is no 
longer lower semicontinuous and we will present counterexamples. 


Example 4.7 (Oscillation patterns). We construct an explicit sequence for which lower semicontinuity 
of the membrane energy ^mem fails. Fix 0 < i? < 1 and Afi = with the parametrization ^ ^ 
Consider a sequence of deformations cpk • ^ defined in polar coordinates of (r, 9) by the 

condition 

(4.29) d^ipkiO = (Rsink^) er{r{(pk{0)^^{^k{0)) + (l - e0{r{(pk{O)^G{^k{O)) ^ 
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where = (cos0, sin0)^, 651 = (—sin0,cos0)^ for given 0^(0). Note that for any k and 9 that 
\de^k{^)\ = 1, so that the transformations are tangentially isometric. We define (pk{0) via two integ¬ 
ration constants tq and 9o for the initialization of r and 0 at ^ = 0. We set 9o = 0 and choose tq such 
that the curve (pk is closed and simple, which imposes tq = r{(pk{0)) = since the first term 

in ( |4.29| ) has zero average. From the second term, taking into account that 651 (r, 9) is independent of r, 
we get the condition 


r27T I I f‘2'Kk I 

27 rro = / (1 — sin^ ^ — j j (l — sin^ C) ^ 

Jo ^ Jo 

where we have applied the change of variables Q = k^. By periodicity the right hand side (an incom¬ 
plete elliptic integral of the second kind with modulus R) is independent of k and thus determines vq. 
The resulting (pk for several values of k are depicted in Figure 

We observe that dop>k{^) in for any 1 < p < 00 (and also weak-* in L^). Therefore, 

the weak VF^’^-limit p) of the p)k is the function defined by p){9) = and obviously not an isometry. 
Assuming 0 < cr < 1 and extending p> along the radial direction to the annulus {I — a < r < 
1 -h cr}, we obtain corresponding deformations given by 

(p{r,e) = ^pk{6) + (r - l)Qi^deipk{0), and </)(r, 0) = ip{e) + (r - l)roe^ = rroe^, 

where (5| stands for clockwise rotation by 7r/2, so that Q^dop)k{9) is the unit outward normal to 

p?k{S^)- Clearly also ^ 0 in on the annulus. We observe that ^mem[0^] = 0, but ^mem[0] > 
0. Hence, ^mem is not weakly lower semicontinuous. 




Figure 3. Explicit oscillations for a simplified model. p)k for R = 0.95, /c = 6, 20,50 


The celebrated Nash-Kuiper theorem |[5niT7]| states that it is possible to uniformly approximate any 
short immersion by isometric ones. Our explicit oscillations around tqS^ is just one example of 
this phenomenon. Notice that a bending term of the type F^bend introduced in our model only compares 
the curvatures of and It therefore does not penalize oscillations, since it does not 

detect the curvature of 0 (AJi) at all. 

Example 4.8. [Lack of rank-one convexity] We present an additional example of a configuration for 
which the integrand of an energy of the type ^mem is not rank-one convex. Rank-one convexity of 
the complete energy density, i.e., , convexity in f G K when composed with the function A^tB for 
any matrix A and any rank one matrix B, is known to be a necessary condition for quasiconvexity 
(I2Q1, Theorem 5.3). Quasiconvexity, in turn, is necessary for weak lower semicontinuity of integral 
functionals in Sobolev spaces (l|20l. Theorem 8.1 and Remark 8.2). 

Let Q = (—2,2)^, and Mi be a closed curve such that Mi H (—1,1) x (0, 2) = (—1,1) x {!}. 
At any point xq G (—1,1) x {!}, the tangential derivatives are just partial derivatives along the first 
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coordinate, yielding 


V„,p(x„) = I)0(w)P(e2) = ( 0 ) ’ 

V,,4>(xofVt,4>{xo) = 0 ^ 

Hence the tangential area distortion measure reduces to 

tr{Vtg 4 ){XQf Vtg(l){XQ)) = det{Vtg(l){XQf Vtg 4 ){XQ) + 62 ® 62 ) 

= {di4>i{xQ))‘^ + {di4>2{xQ))‘^ , 


where 62 = (0,1)^. Defining now the convex function 

F{a, (i) = + dT^ - 2, 

which has a unique minimum with value 0 for a = = 1 , we have that the energy density 

Wf{B) = F (tr(S^S), det(S^S + 62 0 62)) 

has a pointwise minimum, with value zero, whenever Vcj) is such that + ( 5 i 02 )^ = 1- 

Consider now, for 0 < A < 1, the family of matrices 

(431) o)=^(J o)+<'-">(l 

Clearly B{X) is rank one. But we have ^^(^(A)) = + (1 — A)^ + — 2 and therefore 

WpiBm = F{B{1)) = 0, but W^^(5(l/2)) = 
which demonstrates that is not rank-one convex. 


5. Finite element discretization based on adaptive octrees 

We adopt a ‘discretize, then optimize’ approach and consider a finite element approximation and 
optimize for the coefficients of the solution. Since the energy is highly nonlinear and nonconvex, 
we use a cascadic multilevel minimization scheme in which the solution for one grid level is used as the 
initial data for the minimization on the next finer grid. We use adaptive refinement of the underlying 
meshes around the surfaces Af i, Af 2 C (0,1)^ for n = 2,3 (Algorithm[^. 

One of the main characteristics of our functional is the pervasive presence of coefficients depending 
on the deformed position (j){x). Indeed, this is how the functional takes into account the geometry of 
target surface, through the projection P 2 and shape operator 52. From an implementation perspective, 
however, this means that frequently discrete functions have to evaluated at deformed positions. There¬ 
fore, the ability to efficiently search the index of an element containing a given position is of paramount 
importance, so a hierarchical data structure that allows for efficient searching is needed. The model 
only contains first derivatives of the unknown deformation. Hence, multilinear finite elements already 
allow a conforming discretization. For these reasons we use multilinear FEM on octree grids. The 
grids used are such that all of the elements are either squares or cubes of side length h — 2“^, for 
an integer £ to which we refer as grid level of the element. In what follows let us detail the different 
ingredients of the algorithm. 
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Algorithm 1 Cascadic minimization scheme. 

Starting grid: Uniform of level h = 

(j) i — II 

for £ i '^rnin lO '^max dO 

Regenerate di, d 2 on grid by aFMM. 
Compute ni, n 2 , 5i, 52 from di, d 2 . 

(j) ^ L^-CG-descent (0) 

Mark all elements intersecting Mi or Af 2 - 
Refine the grid {h ^ 2“^+^). 
end for 
return (j) 


Multilinear Finite Elements on Octrees. We assume n = 3 for the presentation here. Using an 
adaptive octree grid based on cubic cells leads to hanging nodes (see Figure Q, nodes which are on the 
facet of a cell without being one of its vertices. Enforcing continuity of the finite element functions 
leads to constraints for function values on hanging nodes and these hanging nodes are not degrees of 
freedom. Additionally, to minimize the complexity of the required interpolation rules, the subdivision 
is propagated in such a way that the grid level of neighboring elements sharing a cell facet differs at 
most by one. 

Octrees and the access to degrees of freedom via hashtables. Even though the tree structure gives 
a natural hierarchical structure to the elements of the mesh, maintaining consistent linear indices for 
degrees of freedom, hanging nodes, and elements can be delicate. Consistent rules could be devised to 
maintain consistency with the element octree for a given mesh, but these would not be easy to update 
when the grid is refined. In order to keep track of vertex indices in a simple manner without sacrificing 
efficiency, hash maps (HU, Chapter 11) are maintained to keep track of the indices of degrees of 
freedom, hanging nodes, and cells. The keys used in the hashmap are a combination of a level value I 
and point coordinates as integer multiples ofh — 2“^. These keys uniquely identify nodes or elements, 
with the convention that an element is identified with its lower-left-back corner. Whenever a query for 
a node or cell is made, there are two possible outcomes. If it is already contained in the corresponding 
hash table, a linear index for it can be retrieved. Otherwise, a new entry of the hash table is created 
and the node or cell is given the next unused index. Since we do not require coarsening of the mesh, 
this scheme guarantees a consistent linear set of indices with a computational cost for insertions and 
queries that is, on average, independent of the mesh size. 

Computing distance functions on octrees. In our model we have assumed that the distance functions 
to our surfaces are given. In practice, especially when using adaptive grids, we need to compute signed 
distance functions on such grids. This has been accomplished by a straightforward adaptation of the 
Fast Marching Method on cartesian grids |[57l exploiting the fact that our grids still are subgrids of 
a regular cartesian grid. In the implemented variant hanging nodes are not taken into account for the 
propagation, their values being linearly interpolated to accommodate the constraints needed for con¬ 
formality. The initialization for the distance computation has been performed starting from triangular 
meshes of the surfaces (for n = 3; for n = 2 two-bit segmentation of interior and exterior of the 
curves has been used). The signs of the distance functions have to be computed separately, by detect¬ 
ing which points of the grid are inside (resp. outside) the initial surface data. In our case, they have 
been computed with the provably correct algorithm given in 12. 

Computation of the coefficients. The discretization for the unknown deformation 0, as already men¬ 
tioned, is done by multilinear finite elements. However, the coefficients of our model include first and 
second derivatives of the signed distance functions d^, for the normal vectors and shape operators 
Si (i = 1,2), respectively. The approximations are required to be robust, since they appear in the 
highest order terms of the model. For the normal vectors n^, we compute the projection of the finite 
element derivative of to recover the nodal values of a piecewise multilinear function, followed by a 
orthogonal projection to the unit sphere to restore the constraint |n^| = 1. 
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Figure 4. Hierarchical grids corresponding to the dolphin surfaces (different 2D 
slices in 3D, grid level 8, 178584 DOFs, 1.1% of the amount of DOFs in the full 
grid case) and leaf contours (2D, level 10). 

In the case of the shape operators, our approach is to approximate the distance functions by a 
quadratic polynomial supported on a neighborhood of each point. Given a fixed integer neighborhood 
size r, for each non hanging node (i.e., the neighborhood contains the r closest other degrees 

of freedom of the adaptive grid) the local quadratic polynomial pk is defined as the one minimizing the 
least-squares error 

(Pkixj) - di{xj))^ . 

XjeBr{xk) 

which can be easily computed by inverting a small matrix. The Hessian of at the node x^ is then 
approximated by the Hessian of p^. 

For the computation of matrix square roots and their inverses, we have used the method described 
in ll26l . taking appropriate care to truncate almost-singular matrices, since the resulting square roots 
also appear inverted. 

Minimization strategy. For the minimization at each level, we have opted for a Fletcher-Reeves non¬ 
linear conjugate gradient method ( ll52l . Section 5.2). The gradient of whose computation is 
involved but elementary, was implemented directly. The parameter a is progressively reduced when a 
further feasible descent step is not found, according to an Armijo line search dl52l . Section 3.1). 

6. Numerical RESULTS 

All of our results have been computed on the unit cube = [0,1]^ for the matching of surfaces 
in 3D, and the unit square [0,1]^ for the matching of contour curves in 2D. In practice, we have used 
homogeneous Neumann boundary conditions, since this allows to have relatively large shapes Mi in 
comparison to the size of the domain Vt without creating excessive volume energies near the boundary 
(for the justification we refer to Corollary |4.3| ). However, if the boundary is not fixed, the deformed 
domain ^(fl) is not necessarily contained in Q, so evaluation of coefficients on deformed positions has 
to be appropriately handled numerically. We use a projection of outside position onto the boundary of 
Q for sufficient large dist(A42, dQ). 

For the membrane and the bending energy we use the material parameters A = /i = 1, correspond¬ 
ing to the density 

W{A) = -\Af + -(det^)2 + 

^ 2 4 2 4 
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Figure 5. Behaviour of the optimal (numerical) deformation in the presence of 
strong compression. From left to right: Textured Mi, M 2 , resulting deformed shape 
(/)(Mi) after grid level 7 with our model, and corresponding result after grid level 4 
when P 2 is not present in F^mem- 


In the bending term, the shape operators have been regularized through the truncated absolute value 
function with r = 1. Since we work on the unit cube, this corresponds to a comparatively large 
curvature radius. For the volume term, given that enforcing orientation preservation in a finite element 
framework is a far from straightforward, it is advantageous to work with the simplified version 

Cvol [ W{V^)dx. 

Jn 

We have run the minimization scheme of Algorithm beginning from a uniform grid of level 
^min = 2 or ^niin = 3 with 9^ = 729 nodes, and refined up to ^max = 8 for 3D examples. For 2D 
cases a reasonable range turned out to be ^min = 4,^niax = 10. The finest grids used for two of the 
examples below are depicted in Figure The width of the narrow band was chosen proportional the 
finest resolution of the mesh (a = 2h) since a small value of a clearly produces inaccurate results when 
r]cr is evaluated on coarse grids. However, the constraint = 1 ensures that the overall strength 

of the surface terms E’match, £’mem and F^bend is not affected. The value of the penalty constraint u was 
divided by 10 for each grid refinement, which is justified by Proposition |4.4[ Furthermore, the volume 
weight Cvol was also halved per level to allow for simultaneously higher initial regularization and close 
final matches. Note that this reduction is much slower than that of the matching parameter. 

In all examples, we have used the identity as the initial deformation. It should be noted that although 
the energy is geometric by design, we are using a first-order descent method for its minimization. In 
consequence, an adequate rigid pre-alignment can be beneficial for intricate shapes. Figure shows 
results for the matching of two different dolphin shapes. Our variational approach is highly nonlinear 
and non-convex. Thus, the numerical approximation of the globally optimal deformation depends 
on the initialization of the deformation. Figure shows that the identity deformation as the initial 
deformation is advisable only if the expected optimal deformation is not too large. This is demonstrated 
by applying different rigid body motions to A4i. 

All figures have been produced by deforming the input data (polygonal curve or triangulated sur¬ 
face) via the resulting deformation 0. This is in contrast to deforming the grid and plotting the resulting 
extracted level sets (which effectively visualizes the inverse deformation), as commonly done in the 
registration literature, and also in |[32l . 

Test case. First we present a simple test case to underline the qualitative properties of our model. 
Figure shows a configuration in which a high amount of compression, combined with rotation, is 
required. Our model finds the intuitively correct deformation, but oscillations typical for the lack of 
lower semicontinuity of the underlying energy are induced when P 2 is not used in the membrane and 
bending terms. The bending term assists in matching the curvatures even if the deformation is not 
rigid. Note, however, that for the optimal match the curvature energy E’bend is not expected to vanish, 
as can easily be seen from p.6| ), p.4| ) and the related discussion in Section 
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Figure 6. Detail is added progressively in the cascadic coarse-to-fine scheme. From 
left to right: Textured Af i, resulting deformed shape i) after the computa¬ 

tion on grid level 4,6 and 8, respectively. 


Shape matching applications. We now turn our attention to high resolution examples with real data. 
Figurej^demonstrates the effect of the multilevel descent scheme, in which details are added progress¬ 
ively to avoid spurious local minima. In Figure a high-resolution 2D example is presented. Figures 
and show 3D examples in which the influence of the curvature matching is indispensable to 
obtain shape sensitive matching deformations. For these examples, the shell parameter S was chosen 
quite high, since the curvature matching term E’bend is a major driving force to obtain correct matching 
of geometric features. Table lists the parameter values used, and run times for our implementation. 
We have split the timings between the highest-resolution level and the combined previous ones, since 
in many applications a very high level of detail might not be necessary, thereby significantly reducing 
the required computational effort. 


Fi§- '^min? ^max ^ CyoU ^ at ^min Time, £ < ('^max 1) Time, £ — -^max DOFs at £] 


I 


00 

co' 

0.5 

0.025,0.002 

lh04m 

4h 34m 

695K 

E 


2,8 

0.71 

0.05,0.1 

30m 10s 

lh27m 

313K 



00 

co' 

1 

0.025,0.002 

20m 04s 

50m 50s 

179K 

10 

1 

00 

co' 

0.5 

0.025,0.002 

28m 56s 

lh25m 

408K 


Table 1. Parameters and running times on a workstation with a single Intel Xeon 
E5-1650 CPU (6 cores, 3.2Ghz). Our implementation splits the computation of the 
different terms of the energy and the corresponding derivatives in different threads 
(obtaining a speedup factor 2), but no further parallelization is used. 
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